#----Basics----

rm(list = ls())

# Set working directory
setwd("~/Library/CloudStorage/Dropbox/Whitten_Kagalwala/TS/JOP/Acceptance/KW_Replication/Stationary_UnitRoot")

# Load required libraries
library(tidyverse)
library(haven)
library(stringr)

# Load data
data <- read_dta('stationary_unitroot.dta')
data$time <- c(rep(1, 5), rep(2, 5), rep(3, 5), rep(4, 5), rep(5, 5), rep(6, 5), rep(7, 5))

# Reshape data
phi1_bias <- data %>%
  select(c(time, phi_y, phi_x, starts_with("phi_bias"), starts_with("phi1_bias"))) %>%
  pivot_longer(phi_bias_pa:phi1_bias_adl33, names_to = "model", values_to = "phi1_bias") %>%
  mutate(model = gsub("^ph.*_bias_(.*)", "\\1", model))
phi1_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^phi_(.*)_t1$"), matches("^phi1_(.*)_t1$"))) %>%
  pivot_longer(phi_pa_t1:phi1_adl33_t1, names_to = "model", values_to = "phi1_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
phi1_power <- data %>%
  dplyr::select(c(time,phi_y, phi_x, matches("^phi_(.*)_power$"), matches("^phi1_(.*)_power$"))) %>%
  pivot_longer(phi_pa_power:phi1_adl33_power, names_to = "model", values_to = "phi1_power") %>%
  mutate(model = gsub("^ph.*_(.*)_power", "\\1", model))
phi2_bias <- data %>%
  select(c(time, phi_y, phi_x, starts_with("phi2_bias"))) %>%
  pivot_longer(phi2_bias_adl22:phi2_bias_adl33, names_to = "model", values_to = "phi2_bias") %>%
  mutate(model = gsub("^ph.*_bias_(.*)", "\\1", model))
phi2_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^phi2_(.*)_t1$"))) %>%
  pivot_longer(phi2_adl22_t1:phi2_adl33_t1, names_to = "model", values_to = "phi2_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
phi3_bias <- data %>%
  select(c(time, phi_y, phi_x, starts_with("phi3_bias"))) %>%
  pivot_longer(phi3_bias_adl33, names_to = "model", values_to = "phi3_bias") %>%
  mutate(model = gsub("^ph.*_bias_(.*)", "\\1", model))
phi3_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^phi3_(.*)_t1$"))) %>%
  pivot_longer(phi3_adl33_t1, names_to = "model", values_to = "phi3_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
b1_bias <- data %>%
  select(c(time, phi_y, phi_x, starts_with("b1_bias"))) %>%
  pivot_longer(b1_bias_stat:b1_bias_adl33, names_to = "model", values_to = "b1_bias") %>%
  mutate(model = gsub("^b1_bias_(.*)", "\\1", model))
b1_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^b1_(.*)_t1$"))) %>%
  pivot_longer(b1_stat_t1:b1_adl33_t1, names_to = "model", values_to = "b1_type1") %>%
  mutate(model = gsub("^b1_(.*)_t1", "\\1", model))
b2_bias <- data %>%
  select(c(time, phi_y, phi_x, starts_with("b2_bias"))) %>%
  pivot_longer(b2_bias_ds:b2_bias_adl33, names_to = "model", values_to = "b2_bias") %>%
  mutate(model = gsub("^b2_bias_(.*)", "\\1", model))
b2_type1 <- data %>%
  select(c(time,phi_y, phi_x, matches("^b2_(.*)_t1$"))) %>%
  pivot_longer(b2_ds_t1:b2_adl33_t1, names_to = "model", values_to = "b2_type1") %>%
  mutate(model = gsub("^b2_(.*)_t1", "\\1", model))
b3_bias <- data %>%
  select(c(time, phi_y, phi_x, starts_with("b3_bias"))) %>%
  pivot_longer(b3_bias_adl22:b3_bias_adl33, names_to = "model", values_to = "b3_bias") %>%
  mutate(model = gsub("^b3_bias_(.*)", "\\1", model))
b3_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^b3_(.*)_t1$"))) %>%
  pivot_longer(b3_adl22_t1:b3_adl33_t1, names_to = "model", values_to = "b3_type1") %>%
  mutate(model = gsub("^b3_(.*)_t1", "\\1", model))
b4_bias <- data %>%
  select(c(time, phi_y, phi_x, starts_with("b4_bias"))) %>%
  pivot_longer(b4_bias_adl33, names_to = "model", values_to = "b4_bias") %>%
  mutate(model = gsub("^b4_bias_(.*)", "\\1", model))
b4_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^b4_(.*)_t1$"))) %>%
  pivot_longer(b4_adl33_t1, names_to = "model", values_to = "b4_type1") %>%
  mutate(model = gsub("^b4_(.*)_t1", "\\1", model))
lrm_bias <- data %>%
  select(c(time, phi_y, phi_x, starts_with("lrm_bias"))) %>%
  pivot_longer(lrm_bias_pa:lrm_bias_adl33, names_to = "model", values_to = "lrm_bias") %>%
  mutate(model = gsub("^lrm_bias_(.*)", "\\1", model))
lrm_type1 <- data %>%
  select(c(time, phi_y, phi_x, matches("^lrm_(.*)_t1$"))) %>%
  pivot_longer(lrm_pa_t1:lrm_adl33_t1, names_to = "model", values_to = "lrm_type1") %>%
  mutate(model = gsub("^lrm_(.*)_t1", "\\1", model))
ftest_y_t1 <- data %>%
  dplyr::select(c(time, phi_y, phi_x,matches("^ftest_y_(.*)_t1$"))) %>%
  pivot_longer(ftest_y_adl33_t1:ftest_y_adl22_t1, names_to = "model", values_to = "ftest_y_t1") %>%
  mutate(model = gsub("^ftest.*_(.*)_t1", "\\1", model))
ftest_x_t1 <- data %>%
  dplyr::select(c(time, phi_y, phi_x,matches("^ftest_x_(.*)_t1$"))) %>%
  pivot_longer(ftest_x_adl33_t1:ftest_x_adl22_t1, names_to = "model", values_to = "ftest_x_t1") %>%
  mutate(model = gsub("^ftest.*_(.*)_t1", "\\1", model))
ftest_y_power <- data %>%
  dplyr::select(c(time, phi_y, phi_x,matches("^ftest_y_(.*)_power$"))) %>%
  pivot_longer(ftest_y_adl33_power:ftest_y_adl22_power, names_to = "model", values_to = "ftest_y_power") %>%
  mutate(model = gsub("^ftest.*_(.*)_power", "\\1", model))

#----Ly Bias----
phi1_bias$model <- factor(phi1_bias$model, levels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"), labels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"))
ggplot(phi1_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of y"[t-1])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 8, 3, 4, 5, 6),
                     labels = c("PA", 'Dead Start', "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_phi1_bias.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----Ly Type 1----
phi1_type1$model <- factor(phi1_type1$model, levels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"), labels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"))
ggplot(phi1_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_type1), position = position_dodge(width = 0.65)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 8, 3, 4, 5, 6),
                     labels = c("PA", 'Dead Start', "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_phi1_type1.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----Ly Coverage----
phi1_type1 %>%
  mutate(phi1_coverage = 1 - phi1_type1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_coverage), position = position_dodge(width = 0.65)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 8, 3, 4, 5, 6),
                     labels = c("PA", 'Dead Start', "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_phi1_coverage.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----Ly Power----
phi1_power$model <- factor(phi1_power$model, levels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"), labels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"))
ggplot(phi1_power, aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_power), position = position_dodge(width = 0.65)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Power", x = "T", shape = "", title = expression("Power for the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 8, 3, 4, 5, 6),
                     labels = c("PA", 'Dead Start', "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_phi1_power.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L2y Bias----
phi2_bias$model <- factor(phi2_bias$model, levels = c("adl22", "adl33"), labels = c("adl22", "adl33"))
ggplot(phi2_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi2_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of y"[t-2])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_phi2_bias.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L2y Type 1----
phi2_type1$model <- factor(phi2_type1$model, levels = c("adl22", "adl33"), labels = c("adl22", "adl33"))
ggplot(phi2_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi2_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_phi2_type1.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L2y Coverage----
phi2_type1 %>%
  mutate(phi2_coverage = 1 - phi2_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi2_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability the coefficient of y"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_phi2_coverage.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L3y Bias----
ggplot(phi3_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of y"[t-3])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)")) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_phi3_bias.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L3y Type 1----
ggplot(phi3_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-3])) + ylim(0, 1) + 
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)")) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .)) 
ggsave("Stationary_UnitRoot_phi3_type1.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L3y Coverage----
phi3_type1 %>%
  mutate(phi3_coverage = 1 - phi3_type1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability for the coefficient of y"[t-3])) + ylim(0, 1) + 
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)")) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .)) 
ggsave("Stationary_UnitRoot_phi3_coverage.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----x Bias----
b1_bias$model <- factor(b1_bias$model, levels = c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"), labels =  c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"))
ggplot(b1_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(0, 1, 2, 3, 4, 5, 6),
                     labels = c("Static", "Differences", "PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_b1_bias.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----x Type 1----
b1_type1$model <- factor(b1_type1$model, levels = c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"), labels =  c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"))
ggplot(b1_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(0, 1, 2, 3, 4, 5, 6),
                     labels = c("Static", "Differences", "PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .)) 
ggsave("Stationary_UnitRoot_b1_type1.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----x Coverage----
b1_type1 %>%
  mutate(b1_coverage = 1 - b1_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability for the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(0, 1, 2, 3, 4, 5, 6),
                     labels = c("Static", "Differences", "PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .)) 
ggsave("Stationary_UnitRoot_b1_coverage.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----Lx Bias----
b2_bias$model <- factor(b2_bias$model, levels = c("ds", "gecm", "adl", "adl22", "adl33"), labels =  c("ds", "gecm", "adl", "adl22", "adl33"))
ggplot(b2_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t-1])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(8, 3, 4, 5, 6),
                     labels = c("Dead Start", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .)) 
ggsave("Stationary_UnitRoot_b2_bias.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----Lx Type 1----
b2_type1$model <- factor(b2_type1$model, levels = c("ds", "gecm", "adl", "adl22", "adl33"), labels =  c("ds", "gecm", "adl", "adl22", "adl33"))
ggplot(b2_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(8, 3, 4, 5, 6),
                     labels = c("Dead Start", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_b2_type1.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----Lx Coverage----
b2_type1 %>%
  mutate(b2_coverage = 1 - b2_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability for the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(8, 3, 4, 5, 6),
                     labels = c("Dead Start", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_b2_coverage.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L2x Bias----
ggplot(b3_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t-2])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +   
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_b3_bias.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L2x Type 1----
ggplot(b3_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +    
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_b3_type1.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L2x Coverage----
b3_type1 %>%
  mutate(b3_coverage = 1 - b3_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability for the coefficient of x"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +    
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_b3_coverage.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L3x Bias----
ggplot(b4_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the coefficient of x"[t-3])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)")) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_b4_bias.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L3x Type 1----
ggplot(b4_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-3])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))  +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_b4_type1.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----L3x Coverage----
b4_type1 %>%
  mutate(b4_coverage = 1 - b4_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability for the coefficient of x"[t-3])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))  +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_b4_coverage.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----LRM Bias----
lrm_bias$model <- factor(lrm_bias$model, levels = c("pa", "gecm", "adl", "adl22", "adl33"), labels = c("pa", "gecm", "adl", "adl22", "adl33"))
ggplot(lrm_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Bias", x = "T", shape = "", title = expression("Bias in the LRM of x"[t])) + ylim(-20, 20) + 
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +  
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 3, 4, 5, 6),
                     labels = c("PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_lrm_bias.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----LRM Type 1----
lrm_type1$model <- factor(lrm_type1$model, levels = c("pa", "gecm", "adl", "adl22", "adl33"), labels = c("pa", "gecm", "adl", "adl22", "adl33"))
ggplot(lrm_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x = "T", shape = "", title = expression("Type 1 error rate in the LRM of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +   
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 3, 4, 5, 6),
                     labels = c("PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_lrm_type1.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----LRM Coverage----
lrm_type1 %>%
  mutate(lrm_coverage = 1- lrm_type1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Coverage Probability", x = "T", shape = "", title = expression("Coverage probability for the LRM of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) +   
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 3, 4, 5, 6),
                     labels = c("PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)"))  +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_lrm_coverage.pdf", width = 10, height = 10, units = "in", dpi = 300)

#----F-test, Lags of y Type 1 (sum = true phi_1)----
ftest_y_t1 %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_t1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_ftest_y_t1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----F-test, Lags of y Coverage (sum = true phi_1)----
ftest_y_t1 %>%
  mutate(ftest_y_coverage = 1 - ftest_y_t1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_ftest_y_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----F-test, Lags of y Power (sum = 0 for non-0 phi_1)----
ftest_y_power %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_power), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x= "T", shape = "", title = expression("Power for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_ftest_y_power.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----F-test, Lags of x Type 1 (sum = 0)----
ftest_x_t1 %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_t1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_ftest_x_t1.pdf", width = 25, height = 10, units = "in", dpi = 300)

#----F-test, Lags of x Coverage (sum = 0)----
ftest_x_t1 %>%
  mutate(ftest_x_coverage = 1 - ftest_x_t1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  facet_grid(phi_y~., labeller = label_bquote(phi[y] == .(phi_y), .))
ggsave("Stationary_UnitRoot_ftest_x_coverage.pdf", width = 25, height = 10, units = "in", dpi = 300)
